function [theta,var_eb,invtausq_theta,log_odds,llik_trace] = modQTL_association_prior(Z,M,BURNIN,NEPOCH,theta,invtausq_theta)
%
% [theta,var_theta,invtausq_theta,log_odds,llik] = F(Z,M,BURNIN,NEPOCH,THETA,INVTAUSQ_THETA)
%
% Z              = snp x 1 probability vector
% M              = snp x K GWAS/epigenetic mark matrix
% BURNIN         = burn-in time without sparsity
% NEPOCH         = maximum optimization time
% THETA          = K x 1 weight vector
% INVTAUSQ_THETA = K x 1 sparsity vector
%
% var_theta      = posterior variance of theta
% log_odds       = snp x 1 prior log-odds ratio
% llik           = optimization x 1 log-likelihood trace
%
% code: Yongjin Park, ypp@csail.mit.edu
%

    P_TRUNC = 1e-4;
    TOL = 1e-4;

    Nsnp = size(M,1);
    K = size(M,2);
    assert(Nsnp == size(Z,1));

    % ================================================================
    % initialization
    if nargin > 4,

        theta = THETA;
        invtausq_theta = INVTAUSQ_THETA;

        assert(size(theta,1) == K);
        assert(size(invtausq_theta,1) == K);

    else
        theta = zeros(K,1,'single');
        invtausq_theta = 1e-4*Nsnp*ones(K,1,'single');
    end

    M = single(M);
    M2 = M.^2;

    % ================================================================
    llik_trace = NaN(NEPOCH,1,'single');

    for iter = 1:(BURNIN + NEPOCH),
        % construct quadratic approximation
        eta = M*theta;
        p = arrayfun(@(x) 1/(1+exp(-x)), eta);
        p(p < P_TRUNC) = P_TRUNC;
        p(p > (1 - P_TRUNC)) = 1 - P_TRUNC;
        weight = p.*(1-p); % weight for data points
        response = eta + (Z-p)./weight; % working response

        % take coordinate steps
        for k = 1:K,
            denom_k = sum(weight.*M2(:,k)) + invtausq_theta(k);
            tt = sum(weight.*(response-eta+theta(k)*M(:,k)).*M(:,k))/denom_k;;

            eta = eta + (tt - theta(k))*M(:,k);
            theta(k) = tt;
        end

        var_mle = 1./(M2'*weight);
        var_eb = 1./(1./var_mle + invtausq_theta);

        % estimate shrinkage
        if iter > BURNIN,
            invtausq_theta = 1./(theta.^2 + var_eb);
            invtausq_theta = max(0.1/Nsnp, min(10*Nsnp, invtausq_theta));
        end

        if iter > BURNIN,
            llik = sum(Z.*eta - log(1+exp(eta))) - 0.5*sum(theta.^2 .* invtausq_theta);
            llik_trace(iter-BURNIN) = llik;

            fprintf(2,'NEPOCH = %03d LLIK = %.4e, Vmin = %.4e, Vmax = %.4e\r', ...
                    iter-BURNIN, llik, min(var_eb), max(var_eb));

        else
            fprintf(2,'BURNIN = %03d, Vmin = %.4e, Vmax = %.4e\r', ...
                    iter, min(var_eb), max(var_eb));
        end


        if iter > BURNIN + 5,
            if abs(mean(llik_trace((iter-BURNIN-5):(iter-BURNIN-1)))-llik)/abs(llik+1e-10) < TOL,
                fprintf(2,'Converged\r');
                llik_trace = llik_trace(1:(iter-BURNIN));
                break;
            end
        end
    end

    log_odds = eta;

end
